# ISR Decision-Making ICCTRS paper
# 4/9/2014
# Jonathan Bakdash, jonathan.z.bakdash@us.army.mil

library(data.table)
library(reshape)
library(boot)
library(lsr)

# Read in two data files, one with raw data and the other with the ground truth for asset capabilities on allocation tasks
    ISR_DM_raw_data <- read.csv("M:/My Papers/FY2014/ICCRTS_2014/Data/12.26.13_ISR_DM_raw_data.csv")
    ISR_DM_ground_truth <- read.csv("M:/My Papers/FY2014/ICCRTS_2014/Data/12.26.13_ISR_DM_ground_truth.csv")
  
    #Number of rows 
    rows<-length(ISR_DM_raw_data[,1])
    
    #Create empty table using data.table library
    sdt_table<-data.table(Detect_PredA_Vis = rep(-99, times = rows),  
                          Detect_PredA_IR = rep(-99, times = rows),	
                          Detect_PredA_Radar = rep(-99, times = rows),
                          Detect_Reaper_Vis = rep(-99, times = rows),
                          Reaper_IR = rep(-99, times = rows),
                          Detect_Reaper_Radar = rep(-99, times = rows),
                          Detect_Raven_Vis = rep(-99, times = rows),	
                          Detect_Raven_IR = rep(-99, times = rows),
                          Detect_Global_Vis = rep(-99, times = rows),
                          Detect_Global_IR = rep(-99, times = rows),
                          Detect_Global_Radar = rep(-99, times = rows),	
                          Detect_Shadow_Vis = rep(-99, times = rows),	
                          Detect_Shadow_IR = rep(-99, times = rows)
                          )
    
    #Bind columns for participant and trials to table 
    sdt_table<-cbind(ISR_DM_raw_data[,1],ISR_DM_raw_data[,2],sdt_table)
   
    #ISR_DM_raw_data: Data in columns 3:16, fill in same row, same column on new table
    #Truth table:Data in row 1:8, columns 2:14
      
    sdtvals<-function(rawdata, truth, tabletobefilled) #Classifies responses 
      # Calcuate signal detection values by comparing response value with ground truth value
      # Response | Truth  (0 = absent, 1 = present)
      #     1        1    = Hit (1)
      #     0        1    = Miss (2)
      #     1        0    = False Alarm (3)
      #     0        0    = Correct Rejection (4)   
    {  
      i<-1
      z<-rep(1:8,20) #row value for truth table
      repeat{
        
        {for(j in 1:13) {
          if (rawdata[i,2+j] == truth[z[i],1+j]) #Hit or Correct Rejection
              {
              if(truth[z[i],1+j] == 0) {tabletobefilled[i,2+j] <- "CR"} 
              if(truth[z[i],1+j] == 1) {tabletobefilled[i,2+j] <- "Hit"}
              }      
            
            if(rawdata[i,2+j] != truth[z[i],1+j]) #False Alarm or Miss
              {
              if(truth[z[i],1+j] == 0) {tabletobefilled[i,2+j] <- "FA"}
              if(truth[z[i],1+j] == 1) {tabletobefilled[i,2+j] <- "Miss"}      
              }    
                      }
        }
          i<-i+1
          if(i>rows) {break}
    }
      sdt_tablefilled<<-tabletobefilled #make the filled in table a global var
    }
    
sdtvals(ISR_DM_raw_data, ISR_DM_ground_truth, sdt_table)     

#re-label columns 1 and 2
  colnames(sdt_tablefilled)[1] <-"Participant"
  colnames(sdt_tablefilled)[2] <-"Item"
  #head(sdt_tablefilled)

#Create block variable
#Block = 1: First set of 8 trials (no NIIRS), Block = 2: Second set of 8 trials (NIIRS provided) 
  block<-rep(1:2,each = 8, times = 10) #each = number of trials/items, times = number of unique participants
  sdt_tablefilled<-cbind(block,sdt_tablefilled) #bind block variable to filled table

#Convert data from wide to narrow to make aggregation easier
  data.m <-melt(sdt_tablefilled, id = c("Participant","block", "Item"))
  #head(data.m)


#Aggregated SDT data
#Empty data.table
acc_table<-data.table(Block1 = rep(-99, times = 10),  
                      Block2 = rep(-99, times = 10)
                      )

#Correct rejections + Hits = Accuracy
for(j in 1:10) 
  {
  acc_table[j,1]<-sdtbyparticipantandblock[1,j*2]+sdtbyparticipantandblock[3,j*2]
  acc_table[j,2]<-sdtbyparticipantandblock[1,(j*2)+1]+sdtbyparticipantandblock[3,(j*2)+1]
  }

acc_perc<-data.frame(acc_table/100)


# Robust paired sample t-test
#Function to pass randomly sampled values w/rpl into t.test; on each bootstrap run, extract the t-value, p-values, standard errors, and the effect size 
#Good examples: http://www.statmethods.net/advstats/bootstrapping.html     

#Standard error fct
se<-function(data)
  {
  df<-length(data-1)
  sd<-sqrt(var(data))
  se1<-(sd/sqrt(df))
  return(se1)
  }


ttestboot <-function(data, indices) 
          {
            d <- data[indices,]
            ttest <- t.test(d[,1], d[,2], paired = TRUE) #Passes data into t.test, with boot fct sample is random w/rpl
            tval<-as.numeric(ttest[1]) #Output of t.test[1] = t-value
            pval<-as.numeric(ttest[3]) #Output of t.test[3] = p-value
            se1<-se(d[,1]) #Block 1 SE
            se2<-se(d[,2]) #Block 2 SE
            dval<-cohensD(d[,1],d[,2], method="paired") #Cohen's d for paired design, requires lsr() 
            return(c(tval,pval,dval,se1,se2))
           }
  
  #Bootstrap for t-test
  set.seed(123) #to reproduce results                        
  bootresults<-boot(data = data.frame(acc_perc), statistic = ttestboot, R=1000)
  #t.test(x=acc_perc[,1], y=acc_perc[,2], type=paired)

  #Plot bootstraps
    plot(bootresults, index = 1) #t-values
    plot(bootresults, index = 2) #p-values
    plot(bootresults, index = 3) #Cohen's d
  
  #95% confidence intervals 
    boot.ci(bootresults, index = 3) #Cohen's d
  
#Compare pooled accuracy of block 1 and block 2 to 100% (algo accuracy)
#Take the average of block 1 and block 2 by participant
pooled_acc_table<-c(1:10)
acc_table<-data.frame(acc_table)
  for(j in 1:10) 
    {
    pooled_acc_table[j]<-mean(c(acc_table[j,1],acc_table[j,2]))
    }

pooled_acc_table
mean(pooled_acc_table)

#One-sample t-test: Null = 100 
ttestbootonesample <-function(data, indices) 
  {
  d <- data[indices,]
  ttest <- t.test(data, mu = 100) #Passes data into t.test, with boot fct sample is random w/rpl
  tval<-as.numeric(ttest[1]) #Output of t.test[1] = t-value
  pval<-as.numeric(ttest[3]) #Output of t.test[3] = p-value
  se1<-se(d) #SE
  dval<-cohensD(d, mu = 100) #Cohen's d for paired design, requires lsr() 
  return(c(tval,pval,dval,se1))
  }

#Bootstrap for one-sample t-test
  set.seed(123)                         
  bootresults<-boot(data = data.frame(pooled_acc_table), statistic = ttestbootonesample, R=1000)
  boot.ci(bootresults, index = 3) #d-value



# Exploratory Analyses: Summaries of SDT values
counts <- rep(1, times = 2080)
data.m2 <- cbind(data.m[,1:5],counts)
#head(data.m2)

# By asset
SDTbyAssetNarrow<-aggregate(counts~ value + variable + Participant, data = data.m2, sum)

head(data.m2)

# Means by Participant by item
#Means
SDTbyItemNarrow<-aggregate(counts~ value + Item, data = data.m2, sum)
SDTbyItemPercSummary<-cbind(SDTbyItemNarrow[,1:2],SDTbyItemNarrow[,3]/260)

#SDs: Calculated by participant
SDTbyPsItem<-aggregate(counts~ value + Participant + Item, data = data.m2, sum)
SDTbySD<-aggregate(counts~ value + Item, data = SDTbyPsItem, sd)
SDTbyItemSDSummary<-cbind(SDTbySD[,1:2],SDTbySD[,3]/260)  

#   SDTbyItemNarrow<-aggregate(SDTnumeric~ value + Item, data = data.m2, sum)
#   totals<-aggregate(counts ~ Item, data = SDTbyItemNarrow, sum)
# 
#   SDTPerc<-(SDTbyItemNarrow[,4]/100)
#   SDTbyItemNarrow <- cbind(SDTbyItemNarrow,SDTPerc)
#   head(SDTbyItemNarrow)
#   
#   aggregate(SDTPerc~ value+ Item, data = SDTbyItemNarrow, mean)
#   aggregate(SDTPerc~ value + Item, data = SDTbyItemNarrow, sd)
# 
# Item1<-sum(sdtbyitem[1,2]+sdtbyitem[3,2])/260
# Item2<-sum(sdtbyitem[1,3]+sdtbyitem[3,3])/260
# Item3<-sum(sdtbyitem[1,4]+sdtbyitem[3,4])/260
# Item4<-sum(sdtbyitem[1,5]+sdtbyitem[3,5])/260
# Item5<-sum(sdtbyitem[1,6]+sdtbyitem[3,6])/260
# Item6<-sum(sdtbyitem[1,7]+sdtbyitem[3,7])/260
# Item7<-sum(sdtbyitem[1,8]+sdtbyitem[3,8])/260
# Item8<-sum(sdtbyitem[1,9]+sdtbyitem[3,9])/260
#       
#       PredatorHitCR<-sum(sdtbyasset[1,2:4],sdtbyasset[3,2:4])/sum(sdtbyasset[1:4,2:4])
#       ReaperHitCR<-sum(sdtbyasset[1,5:7],sdtbyasset[3,5:7])/sum(sdtbyasset[1:4,5:7])
#       RavenHitCR<-sum(sdtbyasset[1,8:9],sdtbyasset[3,8:9])/sum(sdtbyasset[1:4,8:9])
#       GlobalHitCR<-sum(sdtbyasset[1,10:12],sdtbyasset[3,10:12])/sum(sdtbyasset[1:4,10:12])
#       ShadowHitCR<-sum(sdtbyasset[1,13:14],sdtbyasset[3,13:14])/sum(sdtbyasset[1:4,13:14])
# 


